pacman::p_load(sf, spNetwork, tmap, tidyverse)In-class Exercise 3
Load Packages
Import Datasets
network <- st_read(dsn="data/geospatial",
layer="Punggol_St")Reading layer `Punggol_St' from data source
`C:\brigittatsai\ISSS626_AY2024-25_T1\In-class_Ex\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
There is Z value in the dataset, drop Z dimension using st_zm() function
childcare <- st_read(dsn="data/geospatial",
layer="Punggol_CC") %>%
st_zm(drop = TRUE,
what = "ZM")Reading layer `Punggol_CC' from data source
`C:\brigittatsai\ISSS626_AY2024-25_T1\In-class_Ex\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
check the output, the Z value has disappeared. The output of the code above interprets the input of the data instead of the output shp file.
Visualising Geospatial Data
Without st_geometry, the plot will be separated:
plot(network)
plot(childcare,add=T,col='red',pch = 19)
pch = point size
add=T, open plot or override the dots to the previous plot
plot(st_geometry(network))
plot(childcare,add=T,col='red',pch = 19)
to put markers using tmap, there are at least 4 ways to do it
tm_symbols(), tm_squares(), tm_bubbles(), tm_dots(), tm_markers()
tmap_mode('view')tmap mode set to interactive viewing
tm_shape(childcare) +
tm_dots(col = "red") +
tm_shape(network) +
tm_lines()tmap_mode('plot')tmap mode set to plotting
leaflet = lightweight mapping package (eg. for mobile app)
Preparing Lixel Objects
length of lixel = 700m -> set to 700 based on study
minimum length of a lixel = 350 -> center point (minimum distance)
lixels <- lixelize_lines(network,
700,
mindist = 350)output: 2645 observation
original (network): 2642 observation
this happens because the last segment is longer than 350
lixels <- lixelize_lines(network,
700,
mindist = 150)if we change minimum distance to 150, the no. of segment will change
samples <- lines_center(lixels)tmap_mode('view')tmap mode set to interactive viewing
tm_shape(lixels) +
tm_lines() +
tm_shape(samples) +
tm_dots(size = 0.01)tmap_mode('plot')tmap mode set to plotting
Performing NKDE
densities <- nkde(network,
events = childcare,
w = rep(1, nrow(childcare)),
samples = samples, # input data
kernel_name = "quartic",
bw = 300,
div = "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)The code chunk below is a way to append the value to sample/ lixel dataframe (concept: left join table, but there is no need for unique identifier):
samples$density <- densities
lixels$density <- densitiesdon’t do sorting (it will change the sequence and the point reference will be inaccurate)
kfunction = accumulative distance
gfunction = ring by ring
The code chunk below results in both kplot and gplot
kfun_childcare <- kfunctions(network, childcare, start = 0, end = 1000, step = 50, width = 50, nsim = 50, resolution = 50, verbose = FALSE, conf_int = 0.05)
To plot only the kfunction, use the following code (change to plotg to plot gfunction):
kfun_childcare$plotk